function [tctab, vtab, ctl] = EPEE_fit(DFILE, CONTROL, sf, steps, stepsel)

global VOLTAGE CURRENT

tctab=[];
vtab=[];
nrec=size(CURRENT);
steps(length(steps)+1)=(DFILE.nr_points-1)*DFILE.rate(1)*DFILE.nr_channel/1000;
disp(steps)
whichstep = stepsel;
lt1=steps(whichstep)/(DFILE.rate(1)*DFILE.nr_channel/1000);
lt2=steps(whichstep+1)/(DFILE.rate(1)*DFILE.nr_channel/1000);
time = make_time(DFILE);
tbase=time(lt1:lt2)-steps(whichstep);
power = 4;
lam=zeros(6,1);
fiterr=zeros(2,1);
vtab=zeros(2,1);
a=zeros(2,1);
t1=zeros(2,1);
t2=zeros(2,1);
t3=zeros(2,1);
% - note that below we see successive fits with the results from the previous fit..
lam=[]; % clear this out first;
k=1;
skipfac = 0; % skip factor (for testing...) 0 implies every trace is done, 1 is every other, etc.
tic;
for j=1:(skipfac+1):nrec(1) % note skip factor...
   i = nrec(1)-j+1;
   s=std(CURRENT(i,1:25));
   m=mean(CURRENT(i,1:25));
   p=max(CURRENT(i,(lt1+5):(lt2-5)));
   initpar = [0,  p, 0.75, 0.5, 20,  1000]; % very initial guess
   vstep = mean(VOLTAGE(i,(lt1+5):(lt2-5)));
   %   disp(sprintf('Entry: %d, peak I = %8.2f, m = %8.2f sd = %8.2f', i, p, m, s))
   if(p > (15*s) & vstep > -80)  % only fit if peak current in the window is well over the noise level, and voltage is positive
      QueMessage(sprintf('Fitting %d (of %d): Vm = %6.1f  Im=%6.1f ', i, nrec(1), vstep, p))
      % reduce the number of points for a first pass fit
      if(j > 1 & ~isempty(lam))
         initpar = lam;
      else % use previous fit for first guess.
         [m, it] = max(CURRENT(i,lt1:lt2));
         tau_rise = tbase(it)*0.63;
         tau_fall = tau_rise*20;
         QueMessage(sprintf('Rise: %8.3f  Fall: %8.3f', tau_rise, tau_fall));
         initpar(3)=tau_rise;
         initpar(4)=tau_fall;
      end
      [lam, fiterr(k)] = EP_fit(tbase, CURRENT(i,lt1:lt2), initpar, power);
      vtab(k) = mean(VOLTAGE(i,lt1:lt2));
      
      a0(k) = lam(1); ai(k)=lam(2); t1(k)=lam(3); t2(k)=lam(4); as(k)=lam(5); ts(k)=lam(6);
      %calculate data to plot
      cur(k,:)=CURRENT(i,lt1:lt2); % cut out only what we need
      am(k)=max(cur(k,:)); % max current in the traces themselves... for comparison
      v=min(tbase):0.25:max(tbase);
      % EPEE 54:
      %fitcur(k,:)=lam(1)+lam(2).*((1-exp(-v/lam(3))).^power).*(lam(4).*(exp(-v/lam(5)))+(1-lam(4)).*exp(-v/lam(6)));
      % EPEE 57:
      % fitcur(k,:)=(lam(1)+lam(2).*((1-exp(-v/lam(3))).^power).*(exp(-v/lam(4)))+lam(5).*(1-exp(-v/lam(6))));
      % EPEE 58:
      fitcur(k,:)=(lam(1)+(lam(2).*((1-exp(-v/lam(3))).^power).*(lam(5).*(exp(-v/lam(4)))+(1-lam(5)))));
      k = k + 1;
   end
end
done=toc;
QueMessage(sprintf('Fitting: %8.3f s', done));

FIT.a0=a0; FIT.ai=ai; FIT.t1=t1; FIT.t2=t2; FIT.as = as; FIT.ts = ts;
FIT.err = fiterr; FIT.v=vtab; % FIT.func=sprintf('a0+ai*(1-exp(-t/t1))^%d*as*(exp(-t/t2)+(1-as))', power);

CONTROL(sf).FIT=FIT; % add fit to results.
ctl = CONTROL;

h = findobj('Tag', 'EPEE_FIT'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'EPEE_FIT', 'Name', 'EPEE Time couse fits', 'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...

orient landscape
% set(gcf, 'FontName', 'Arial');

fsize = 7;
msize = 3;
subplot('position', [0.15 0.32 0.30 0.20]);
plot(vtab, ai, '-ko', vtab, am, '-k', vtab, a0, '-bs', vtab, as, '-r^', 'MarkerSize', msize);
set(gca,'FontSize', fsize, 'FontName','Arial' );
%xlabel('V (mV)', 'Fontsize', 10);
ylabel('A (pA)', 'FontSize', fsize, 'FontName', 'Arial');

subplot('position', [0.15 0.075  0.30 0.20]);
plot(vtab, t1, '-ko', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
xlabel('V (mV)', 'FontSize', fsize);
ylabel('{\tau}{_1} (ms)', 'FontSize', fsize);

subplot('position', [0.55 0.075  0.30 0.20]);
plot(vtab, t2, '-gx', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
xlabel('V (mV)', 'FontSize', fsize);
ylabel('{\tau}{_2} (ms)', 'FontSize', fsize);

subplot('position', [0.55 0.32  0.30 0.20]);
plot(vtab, fiterr, '-gv', 'MarkerSize', msize);
set(gca, 'FontSize', fsize);
%xlabel('V (mV)', 'FontSize', 10);
ylabel('Error (pA)', 'FontSize', fsize);

subplot('position', [0.15, 0.57, 0.7, 0.35]);
set(gca, 'FontSize', fsize);
%set(gca, 'Xscale', 'log');
line(tbase, cur,...
   'color', 'black',...
   'Marker', '.',...
   'MarkerSize', 1,...
   'LineStyle', 'none');
line(v, fitcur,...
   'color', 'red');

title(sprintf('File: %s  Recs: [%d-%d] Window: %7.1f-%7.1f\nEPEE fits', ...
   DFILE.filename, DFILE.frec, DFILE.lrec,  steps(whichstep), steps(whichstep+1)) ,'FontSize',7);

return

function [time] = make_time(DFILE)
time=0:(DFILE.rate(1)*DFILE.nr_channel)/1000:((DFILE.nr_points-1)*DFILE.nr_channel*DFILE.rate(1))/1000;
return



function [lam, fiterr] = EP_fit(x_data, y_data, initpar, alpha)
% A current activation and inactivation time course fitting
% model=57;	%Modified EPEE ( m^k*h + n) lam(4) is n amplitude; lam(6) is n activation time constant
model=58;	%Modified EPEE ( m^k*(h*(1-p))) lam(5) is fraction of inactivating channels

%x_pass1 = [min(x_data):0.25:20,25:5:max(x_data)];
%y_pass1 = interp1(x_data, y_data, x_pass1);
%%FitData(:,1)= x_pass1';
%FitData(:,2)= y_pass1';
FitData(:,1) = x_data';
FitData(:,2) = y_data';
[m, i] = max(y_data);
tau_rise = x_data(i)/3;
tau_fall = tau_rise*30;
%tau_fall = exp_taylor(x_data, y_data);
%fit transient current
imax=max(FitData(:,2));
imax2 = 2*imax;
% pars:    Y0      A1       T1        t2    an      Tn
%initpar = [0,    imax*0.5,    tau_rise,    tau_fall,   0.5,    100]; % very initial guess
pmask =   [1,       1,      1,    1,   1,      0];
lbound =  [-imax,   0,    0.1,    0,   0,      1];
ubound =  [imax, imax2,    20,  100,   1,    2000];
order=length(initpar)-1;
maxiter = 1500;
beta = 0;
lam = initpar;
watchon;
%[fpar2, chisq2, niter2, volt_fit2(i,:)] = mrqfit('exponential', [amp0 amp1/2 tau1 amp1/2 tau1/8], time_fit, vsmo(rec, t_list(rec,:)), [], ...
%			[], [], [], 50, []);

[c,lam]=curve_fitting(FitData(:,1), FitData(:,2), 'levenberg','cubic',model,order,initpar, pmask, lbound, ubound, alpha, beta, maxiter);

% now repeat, using the initial values to seed the fit of the full data set for another 100 iterations
%initpar = lam;
%clear FitData;
%FitData(:,1)= x_data';
%FitData(:,2)= y_data';
%maxiter=250;
%[c,lam]=curve_fitting(FitData(:,1), FitData(:,2), 'levenberg','cubic',model,order,initpar, pmask, lbound, ubound, alpha, beta, maxiter);

watchoff;
%now plot it
t = FitData(:,1); 
z = FitData(:,2);
f=fit_func(lam, x_data, y_data, model, alpha, 0);

result = sprintf('Y0=%7.1f A=%7.1f T1=%7.3f T2=%7.3f T3=%7.3f F=%5.3f err=%4.3f',...
   lam(1), lam(2), lam(3), lam(5), lam(6), lam(4), norm(f));
fiterr=norm(f);
%disp(result)

return;


function [tau, A] = exp_taylor(x, y);
% find exp tau based on taylor series expansion and polynomial fit to falling part of trace.
% simple, yet effective. Works for only one exponential..
%

[m, i1] = max(y); % find max point
i2=length(y);

p = polyfit(x(i1:i2),y(i1:i2),9); %fit with 7th order polynomial
A = p(8); % get amplitude here
tau = -p(8)/p(7);
QueMessage(sprintf('ExpTaylor: A=%7.3f tau=%7.3f', A, tau));
return;

